Steady shear flow thermodynamics based on a canonical distribution approach 
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A non-equilibrium steady state thermodynamics to describe shear flows is developed using a 
canonical distribution approach. We construct a canonical distribution for shear flow based on the 
energy in the moving frame using the Lagrangian formalism of the classical mechanics. From this 
distribution we derive the Evans-EIanley shear flow thermodynamics, which is characterized by the 
first law of thermodynamics d£ = TdS — Qd-y relating infinitesimal changes in energy £, entropy S 
and shear rate 7 with kinetic temperature T. Our central result is that the coefficient Q is given by 
Helfand's moment for viscosity. This approach leads to thermodynamic stability conditions for shear 
flow, one of which is equivalent to the positivity of the correlation function of Q. We emphasize 
the role of the external work required to sustain the steady shear flow in this approach, and show 
theoretically that the ensemble average of its power W must be non-negative. A non-equilibrium 
entropy, increasing in time, is introduced, so that the amount of heat based on this entropy is equal 
to the average of W. Numerical results from non-equilibrium molecular dynamics simulation of 
two-dimensional many-particle systems with soft-core interactions are presented which support our 
interpretation. 

Pacs numbers: 05.70.Ln, 83. 50. Ax, 05.20.Jj, 83.60.Rs 



I. INTRODUCTION 

The great success of thermodynamics as a physical the- 
ory to describe various equihbrium phenomena has stim- 
ulated attempts to generalize it to a theory applicable 
to macroscopic time-dependent phenomena, namely to 
a non-equilibrium thermodynamics. Many efforts have 
been devoted to this subject, and led to some propos- 
als for non-equilibrium thermodynamics, for example, 
the classical irreversible thermodynamics the ratio- 
nal thermodynamics 0, Q and the extended irreversible 
thermodynamics 0|. Recently, a non-equilibrium ther- 
modynamics, which tries to give more rigorous predic- 
tions by restricting its applied field into non-equilibrium 
steady states, is also discussed 7]. 

Shear flow is a typical example of non-equilibrium 
steady phenomena. For a constant velocity gradient it 
has a steady current (the shear stress), and has many 
applications in the investigation of rheological properties 
of materials Such models have been widely used 

to calculate the shear viscosity, whose shear rate depen- 
dence is still actively discussed ^> ^| • The ap- 
parent existence of critical phenomenon, appearing as a 
transition from a uniform bulk phase to an organized 
string- like phase, is shown at high shear rate [l^ IirL [lfil | . 
It is used as an application of non-Hamiltonian molec- 
ular dynamics called the non-equilibrium thcrmostated 
dynamics |17| . by which some non-equilibrium dynam- 
ical properties such as the conjugate pairing rule for 
the Lyapunov spectrum [T^ |3 HO; 113 
tuation theorem jU were discussed. Shear flow has 
also been described by the Bhatnagar-Gross-Krook ki- 
netic equation, which is a simplification of the Boltzmann 
equation, and using this equation, transport coefficients 
and hydrodynamic modes were calculated HE Un- 



steady shear flow is a spatially homogenous and time- 
independent phenomenon, so a simple description com- 
pared to general non-equilibrium phenomena may be ex- 
pected. However, although it is such a simple steady non- 
equilibrium phenomenon, a convincing thermodynamic 
description of shear flow is not known yet. 

A proposed non-equilibrium steady state thermody- 
namics to describe shear flow was given by Evans and 
Hanley (27k. ,22,, .2£, . It expresses the first law of ther- 
modynamics for the shear flow by adding the term 
expressing the response to a shear rate 7, namely 



d£ = TdS + ^d-f 



(1) 



as a relation among infinitesimal quasi-static changes 
of internal energy £, entropy S and the shear rate 7 
with temperature T and the coefficient ^ defined by 
^ = dE/d^\s [13. As a conceptual feature, the Evans- 
Hanley thermodynamics is characterized by the fact that 
the shear rate is an external parameter chosen as an 
additional variable to describe non-equilibrium effects. 
This is analogous to the choice of variables in equilib- 
rium thermodynamics where the variables are chosen 
as parameters manipulated externally, for example, the 
temperature and volume, etc. This choice of thermo- 
dynamic variables has the advantage that observables 
in this approach are rather easy to access by experi- 
ments and computer simulations. This feature also dis- 
tinguishes the Evans-Hanley thermodynamics from some 
other non-equilibrium thermodynamics, in which local 
quantities related to conserved quantities, such as the 
local momentum density, are chosen as additional vari- 
ables to describe non-equilibrium effects, because they 
change slowly with time and are consistent with the phe- 
nomenological equations of hydrodynamics. Compared 
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with such a general formahsm of non-equihbrium ther- 
modynamics, the Evans-Hanley thermodynamics gives a 
much simpler description, as its applied field is restricted 
to steady shear flow. On the other hand, one of the 
problems in the Evans-Hanley shear flow thermodynam- 
ics was that a clear physical meaning for the coefhcient 
f , especially its microscopic expression, was not known, 
so that one has not had clear experimental or numerical 
evidence to support this thermodynamics. On this point, 
Ref. tried to calculate numerically a value of ^ by 
introducing a non-equilibrium entropy in a low density 
system. Also recently, Ref. 33j discussed a phenomeno- 
logical expression for using the Maxwell model for shear 
flow. 

Another important aspect of non-equilibrium thermo- 
dynamics is its construction from a solid statistical me- 
chanical foundation. Some attempts in this direction 
have been discussed usin g th e non-equilibrium canoni- 
cal distribution approach IsM ISa IsS IsS ^dl ^tnd the 
projection operator approach |40l l4ll l42l l43l lH . l45l| and 
so on. As one such approach, the non-equilibrium canon- 
ical distribution approach justifies the response formula 
for thermodynamic perturbations, which is different from 
a mechanical perturbation expressed as a change in an 
external parameter appearing explicitly in the Hamilto- 
nian It uses, in principle, the distribution /(F) of 
canonical type 



/(r)=S-icxp<^-/3 



H{r) 



;A*c(rMo(r) 



(2) 



where S is a normalization constant, (3 is the inverse tem- 
perature, H(T) is the Hamiltonian as a function of the 
phase space vector F, and ^a(r) and Aa{T) are pairs 
of conjugate variables whose forms depend on the non- 
equilibrium phenomena under consideration. In many 
cases, the distribution /(F), and therefore the functions 
^q(F) and Aa{T), a = 1,2,- •• ,n, is introduced based 
on the "local equilibrium assumption" [s^, although it 
is not always necessary. Using the distribution /(F) and 
the Liouville operator L, we calculate average quantities 
as the ensemble average under time-evolution of the dis- 
tribution function 



/(r,t) = exp{-zL(t-to)}/(F), 



(3) 



which may be regarded as the distribution function at 
time t evolved from the initial canonical distribution 
function /(F) at time to. In this way, we can derive 
the thermal response formula for viscosity, thermal con- 
ductivity and so on [s^, IH, ■ This approach was 
generalized to non-Hamiltonian systems such as the SUod 
equation for shear flow systems with isokinetic thermo- 
stat ^3, ; and was used to calculate some thermal 
quantities, such as the specific heat of non-equilibrium 
steady states |H3 1^ ■ However it is still an open prob- 
lem to construct the Evans-Hanley shear flow thermody- 
namics from this non-equilibrium canonical distribution 



approach based on distributions of this type 10). Usually, 
the thermodynamic relation for the first law of thermo- 
dynamics in the canonical distribution approach is in- 
troduced based upon the local equilibrium assumption, 
but the form as the first law of thermodynamics for 
shear flow cannot be justified by this assumption, be- 
cause Eq. including the term ^c?7 expressing the non- 
equilibrium effect of the shear flow cannot be attributed 
to an equilibrium thermodynamical relation even if we 
consider a very small portion of system. 

The principal aim of this paper is to derive the Evans- 
Hanley shear flow thermodynamics based on the canon- 
ical distribution approach to non-equilibrium steady 
states. First we show that the canonical distribution for 
shear flow is represented by the distribution Q in the 
case where n — 1, /iQ(F) — 7 and Aa{T) — — Q(F). Here 
(5(F) is the Helfand moment of viscosity and its time- 
derivative is connected with the off-diagonal component 
of the pressure tensor. The canonical distribution given 
here is different from the canonical distribution based 
on the local equilibrium assumption, because the "non- 
equilibrium term" —jQ{T) in the distribution /(F) can- 
not be neglected even if we take the system size to be 
very small. So to derive the canonical distribution for 
shear flow, we introduce the Hamiltonian for the mov- 
ing frame which follows the steady global current, then 
using the Lagrangian techniques of classical mechanics, 
and the fact that the quantity H{T) + J2a^i Aia(F)Aa(F) 
in the distribution ^ should correspond to this moving 
frame Hamiltonian. This procedure, which is one of the 
important points of this paper, gives a systematic way to 
choose the function /^Q(r)AQ(F) for the canonical 

distribution approach to non-equilibrium steady states. 
To justify this procedure concretely, we show that by ap- 
plying it to the rotating system we can obtain the well 
known canonical distribution for the rotating system, and 
therefore its thermodynamics. 

Second we show that our approach is consistent with 
the linear response formula for viscosity, which can be 
derived from the ensemble average of the shear stress 
using the distribution Q. To interpret the derivation of 
this response formula we emphasize the role of the work 
required to sustain the steady shear flow. We introduce 
a non-equilibrium entropy, which increases in time, and 
show that the heat based on this entropy has the same 
magnitude as the power needed to sustain the shear flow. 

Third we derive the form Q of the first law of thermo- 
dynamics for shear flow from our canonical distribution 
approach, and show that the quantity ^ in the form 
is given by ^ — —Q with Q being the ensemble aver- 
age of the Helfand moment for viscosity Q(F). We also 
discuss the thermodynamic stability condition for shear 
flow, which leads to the positivity of the correlation func- 
tion of Q(F) as well as the positivity of the specific heat. 

Finally we present some numerical calculations of 
many-particle systems with soft-core interactions to sup- 
port our thermodynamic interpretation of steady shear 
flow. Here, we use the Sllod equations with the isoki- 
netic thermostat and Lee-Edwards boundary condition 
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■ In these simulations we show the shear rate depen- 
dences of the average of Helfand's moment of viscosity, 
its correlation function and the work needed to sustain 
the shear flow. 



this frame change — > j:(mov) ^ Inserting the veloc- 
ity v^*"'^-' — vj"'"^-* +V{q^j) into Eq. the Lagrangian 
l(niov) ^ ^(mo«)^^(mo^)^q) ^^le systcm in the moving 
frame J^t^o") is given by 



II. CANONICAL DISTRIBUTION FOR 
STEADY FLOWS 

A. Moving Frame and Energy 

Systems discussed in this paper are steady flows, in 
which the global velocity distribution of the flow is given 
by a time-independent function V(r) at the position r in 
the inertial frame We consider a system which 

consists of N particles with the same mass m and is de- 
scribed by classical mechanics (without magnetic fields) . 

In this section we use the Lagrangian formalism of clas- 
sical mechanics to compare quantities in different frames. 
The Lagrangian formalism is a direct consequence of the 
frame-independent principle of least action S J^^ dtL = 
for fixed values of positions at times ti and t2, where L 
is the Lagrangian as a function of particle positions and 
velocities. 

In the inertial frame jf(™'=) the Lagrangian L^™'^') = 
X(»"'=)(v(™'^\ q) is given by 



Jine) 



f/(q) 



(4) 



1 ^ , 

J- \ (mov) 



V(q,) -t/(q) 



(8) 



r r (rnov'\ ( (rnov) (rnov) (mov)\ 

as a function of v*-™°^^ = (v^ ••• ,v]y ') 

and q. Equation © leads to the momentum p(™°^) = 



as 



3 

{ine) 



v(q,; 



{ine) 



(9) 



Therefore the momentum is independent of choice of the 
frames J^t™"") and and hereafter we use the no- 



tation p = (pi, p2 • • • , Pat) = p 



{mov) 



momentum, and also use the notation T = (p, q) for the 
phase space vector. On the other hand, the Hamiltonian 
j^(mov) ^ ^^mov) -^^ ^j^g moving frame is given 

by 



as a function of v^™''-' = (vj^™*^ , ^r'. 



(ine) (ine) 



q = (qi, q2, • • • , qw j, where v^*"'^'' and q^ are the veloc- 
ity and the spatial position of the j-th particle, respec- 
tively, and U{q) is the potential function as a function 
of q. Using the definitions p^™'^) = 9L(™'=)/5v(™'^) and 



) and 



^(««e)(p(.me)^q) EE p(»«e) . v^™'^) - L^™'^) , the Lagrangian 



Q leads to expressions for the momentum p^^' 
Hamiltonian ij'* 



as 



and the 



(5) 



N 



(10) 



using the definition i/('"°^)(r) = p • v^""") - L^™"") and 
Eq. ||HJl. It is essential to note that although the momenta 
pimov) g^^^ pizne) ^j.g ^^^^^^ ^j^g Hamiltonian iJ (™°^) (r) 

in the moving frame jf(™°'") is different from the Hamil- 
tonian fj(tne)^Y) in the inertial frame and their 
difference is proportional to the global current distribu- 
tion V(qj) ,53<| . 



ff("")(p('""\q) 



2m 



(me) 



C/(q). 



(6) 



Here we note that the Hamiltonian is a function 

of p(»e) and q, namely = ij(™'=)(p("'=), q). 

Now, we introduce the velocity v^™""^ of the j-th par- 
ticle in the moving frame J^i"^""")^ which is connected to 
the velocity vj™'^'' in the inertial frame jf(™^) by 



(mov) (ine) 



The quantity v^™""'' is often referred to as the thermal 
velocity of particle j. The position q is invariant under 



-V(q,) 



(7) 



B. Canonical Distribution 

The central assumption of this paper is that in the 
moving frame defined by the global velocity dis- 

tribution V(r) the system can be regarded as an equilib- 
rium state. It is important to note that this assumption 
is not obvious, because generally a global flow causes 
some local effects such as string phases in shear flows 
[T^ lisl [T^ [r^ and possibly turbulent phases and so on, 
which can destroy this assumption. However these phases 
generally occur in far-from equilibrium states so we may 
expect that our assumption is satisfied in a regime near 
equilibrium. Under this assumption, we introduce the 
canonical distribution 
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/(r) = s-icxp{-/?i7(™°'')(r)} 

= exp { -/? 



(11) 



N 



(12) 



for steady flow, where /3 is the inverse temperature 
l/{kBT) with the Boltzmann constant fcs and the 
temperature T, and S is the partition function S = 
/ dr exp {— /37?''"°'"^(r)}. For simplicity we use units 
so that ks — ^ hereafter in this paper. It should 
be noted that the function f{T) is the distribution of 
p(me) q ^j^g incrtial frame jf(™'=)^ as well as 

the distribution of p^™"") and q in the moving frame 
jr(mou) ^ because of the identity of the momenta implies 
r = (pt^-^^q) = (pf^^^^q). 

The canonical distribution (|12|l is different from that 
obtained from the " local equilibrium assumption" , which 
is popular in many texts on non-equilibrium statistical 
mechanics. In this approach the canonical distribution 
function g{T) is chosen as 



= ^ cxp 



exp 




E^|p-"^v(q^)l' + f^(q 

i=i 

^ 1 

if("-)(r)+^ ™|V(q 




with S = /drexp{-/3Ef^Jp - ,7iV(q^.)lV(2m) + 
t/(q)]}. This type of distribution was actu ally used to 
calculate the shear viscosity [M IMIlHHIlil , and its 
localized version was used as a canonical distribution un- 
der the local equilibrium assumption (See for example, 
Ref. js^ ) • However it is not consistent with the thermo- 
dynamics of rotating systems discussed in the next sub- 
section, so the difference between the two distributions 
(|12|1 and H13|l is crucial to the subject of this paper, that 
is a statistical foundation for steady state thermodynam- 
ics. Another problem with the distribution .g(r) is that 
it does not take into account the incrtial force. Despite 
these problems, one may also notice that the deviation of 
the distribution H13|) from the distribution H12|) is of or- 
der 0(|Vp) in the global velocity, and near equilibrium 
it may give small effects compared with the effects given 
by the second term on the right-hand side of Eq. IIUI) 
which is of order OdVI). 



III. ROTATING SYSTEMS AND THEIR 
THERMODYNAMICS 

Before discussing the canonical distribution approach 
to shear flow systems, in this section we give a derivation 
of the well known canonical distribution and the thermo- 
dynamics for uniformly rotating systems, based on the 
formalism given in Sec. ^ The rotating system is a 
typical example to which our approach can be applied, 
and a comparison between the rotating system and the 
shear flow system gives a useful insight into the canonical 
distribution approach and the thermodynamics of shear 
flow systems which will be developed in Sees. II VI and FVl 



A. Canonical Distribution for Rotating Systems 

We consider a rotating system with a constant angu- 
lar velocity vector u). We assume that the Hamiltonian 
i?^™^^(r) is invariant under rotation about the axis of 
rotation. In this section the origin of the spatial coordi- 
nates and the axis of rotation is taken at the center of 
mass of the system. Under these conditions the global 
velocity distribution function V is given by 



V(qj) =a; X qj 



(14) 



where x is the usual vector product. This global velocity 
distribution can be sustained without any external effect 
in isolated systems, because of the conservation of the 
total angular momentum. 

Using Eq. (fT^ . relation lfTU|l is rewritten as 



(15) 



with the total angular momentum M(r) = X^jLi 1j ^ Pj- 
In comparison with shear flow systems discussed in Sec. 
IIVI it is important to note that the total angular mo- 
mentum M(r) is conserved in the inertial frame J^("^<^) ^ 
namely 



jL(™'=)M(r) = 0, 



(16) 



where L^™*^) is the Liouville operator defined by iL(*"^)X 

^ (aij(»e)(r)/ap) . {dx/dq) - (ai?(™'=)(r)/aq) • 

(dX/dp) for any function X of T. Using Eqs. l(T^ and 
()15|l the canonical distribution for the rotating system is 
represented as 

/(r) = exp [ijt"''^' (r) - a; • M(r)] } . (17) 

The distribution H17(l is stationary, namely iL^'^'''^'^^ f iT) = 
ii:("°")/(r) = 0, in both frames and 
Here L^™"") is the Liouville operator in the moving frame 
j:{mov) ig defined by iL^'"™)^ = (aij(™°")(r)/5p) • 
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(dX/dq) ~ (ai? (r) /9q) • (dX/dp) for any functfon 
X{T), similarly to The distribution (|T7J| has the 

general form l^Jl of the canonical distribution in the case 
that H{T) = H^''"^^T), h = J, and Aj{r) is the compo- 
nent of M(r), and /ij(r) is the component of u) in the 
d-dimensional system. 

It is valuable to note that from the canonical dis- 
tribution (|17|l we can derive the distribution function 
y(v('"°"),q) for the position q and the velocity v*^™""^ 
in the moving frame jf(""'") as 




(rnov) 



+ U{q) + u{q) 



(18) 

using Eq. where u(q) is given by u(q) = 

— 5™ I'^l^ ^'j with Vj = |a; X qj|/|a;|. The function 

j'(v(™°"), q) is the distribution for a rotating system in- 
cluding explicitly the effect of the centrifugal potential 
w(q). We cannot derive the distribution p8|l from the 
distribution 



B. Thermodynamics for Rotating Systems 

In this paper we use the notation X for the ensemble 
average of any function X(T) using the canonical distri- 
bution /(r), namely 

X = J dTX{T)fiT). (19) 
Using this notation and Eq. H15() we obtain the relation 



(20) 



which connects the average energies ijC™"") and in 
the two different frames J^i"^"") and ^ respectively. 

Now we introduce the observable S(T) corresponding 
to entropy as 



p{mov) ^ jj{mov) _ g 

= -TlnS 



(24) 
(25) 



where we used Eq. (|22|l to derive Eq. (|25|l . Similarly 
the free energy F*^™"^) in the inertial frame 
introduced as 



p{ine) ^ f{{ine) _ pg 

= -T I InS - 



dlnE 
du3 



(26) 
(27) 

(28) 



using Eq. (|^ . Therefore the free energies 

can be calculated directly from the partition func- 
tion S. 

Noting the definition of the partition function S = 
/ dT ea;p{-/3i/("™)} = j dV exp{-P - uj ■ M]} 

including the two parameters T{= f3~^) and to explicitly, 
Eq. H25|) implies that the free energy pi"^°") is a func- 
tion of the temperature T and the angular velocity lu: 
pimov) ^ f'-"''"'\T,ui), and we have 



dp 



dlnE 
d(3 



91nS 



Eqs. and l|5n|l are summarized as 



= -/3M 



ppimov) ^ H(""'^)d(3- (3M-dLJ, 



(29) 



(30) 



(31) 



where we have used the notation that dX is the infinites- 
imal change in the quantity X. Inserting Eqs. H24() and 
d(3 = -T-'^dT into Eq. (EB we obtain 



^pimov) ^ _g^p -M-duj. 

Noting Eq. (|27|l . we also obtain 



(32) 



5(r)^-in{/(r)}, 

so that the entropy S is given by 



S = ln5 + /3iJ(" 



(21) 



(22) 
(23) 



where we used Eq. p7|) . The free energy F^™'°'"'> in the 
moving frame J^i"^"") is also introduced as 



^pizne) ^ ^SdT + U}- dM. 

Equation (|32() is also equivalent to 



= TdS -M-duJ, 



^^(,ne) = TdS + LJ ■ dM, 



(33) 



(34) 



(35) 



noting the relations (|20|l and 124|) . The second term on 
the right-hand side of Eq. (|34() can also be derived from 
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the relation dH^'"'°''^T)/du; = -M{T) derived from Eq. 
(Unj, therefore SiJC""") /du:\-g = -M under an adiabatic 
process. From Eq. (|34|l . the energy 7J('"o'") in the moving 
frame J^^"^°^) is regarded as a function of S and u;, while 
the energy in the inertial frame ^(*"^) is a function 

of 'S and M by Eq. The relation ^ is the first 

law of thermodynamics for the rotating system, which is 
well known |5^ . 

Using Eq. jSSJ we obtain dS/dH^]^ = l/T and 
dS / dyVq^i—^ — u>/T hy regarding the entropy 5 as a 

function of and M. We notice that the thermo- 

dynamic variable conjugate to the averaged angular mo- 
mentum M is the inverse temperature times the angular 
velocity oj/T, like the fact that the thermodynamic vari- 
able conjugate to the energy TJf™*^) is the inverse tem- 
perature 1/r. 

After all, thermodynamic functions such as the free 
energy 

pimov) 

are calculated from the partition function 
S, and by combining them with the first law of thermo- 
dynamics we can calculate thermodynamic quantities, for 
example, the total momentum M = — 9F(™°^)/9a;|T and 
the entropy S = —dF^™-°^^/dT\ijj, and their relations in- 
cluding the equation of state, are as in equilibrium ther- 
modynamics. 



CANONICAL DISTRIBUTION APPROACH 
TO SHEAR FLOWS 



IV. 



In Sec. mil we discussed the construction of the well- 
known thermodynamics of rotating systems from the 
canonical distribution H12|l . Clearly we can carry out a 
similar procedure for shear flow systems, as will be shown 
in Sec. However it is not clear that the thermody- 
namic relations obtained by such a procedure correspond 
to the shear flow thermodynamics proposed by Evans 
and Hanley. This section is devoted to discussing this 
point. We compare the shear flow system with the rotat- 
ing system discussed in Sec. IIIII and emphasize not only 
similarities but also differences between these two sys- 
tems. This difference leads to the necessity to consider 
the work required to sustain steady currents, and plays 
an essential role in deriving the linear response formula 
for viscosity in the shear flow system. 



A. Shear Flows and Helfand's Moment of Viscosity 

We consider the shear flow system in which the global 
current distribution is given by 



(36) 



where qjy is the y-component of the spatial coordinate 
Qj of the j-th particle, and '\x is the normalized vector 
pointing to the positive a;-direction. Here, 7 is the shear 
rate, namely the constant gradient of the x component 



of the local velocity as a function of y, and is assumed to 
be a position- independent constant. 

The shear flow system was proposed to describe fluid 
filled between the two plates which move at different 
speeds, and it is frequently used to calculate viscosity. 
The viscosity is calculated as the linear coefficient of the 
average of the a:y-element P^yiT) of the pressure tensor 
Pap{T) defined by 



N 



9C/(q) 



(37) 



as a function of the shear rate 7. Here V is the vol- 
ume of the system, pja (qja) is the a-th component of 
the momentum pj (the spatial coordinate q^) of the j- 
th particle. If particle-particle interactions in the sys- 
tem are given by a two-body interaction only, namely 
the potential C/(q) is expressed in the form C/(q) — 
(1/2) (/)(|qj — qfel), then Eq. f^7\i can be rewrit- 

ten as 



^0/3 (r) 



1 ^ 



1 



1 ^ 



k=l 



(38) 

with Tjfca = ~{l/2)d(j){\qj — qjt\)/dqja interpreted as the 
a-th component of the force acting to the j-th particle 
by the A:-th particle. The pressure tensor Paf}{T) comes 
from the balance equation for the momentum 17]. 

For the case of the global velocity distribution H36|l . 
Eq. (fTn|l is given by 



where Q{r) is defined by 



Q(r) 



N 



(39) 



(40) 



It is essential to note that the quantity Q{T) is connected 
to the shear stress P.j;y{r) as 



VP.,(r), 



(41) 



narnely the quantity Q{T) is Helfand's moment of viscos- 
ity |53,|5a|- (In the references, the name "Helfand's mo- 
ment of viscosity" is used for the quantity {TV)'^^'^Q{T), 
but in this paper, for convenience we use this name for the 
quantity Q itself without the factor {TV)~^^^.) Helfand's 
moment of viscosity is used to calculate the viscosity by 
analogyjvith the Einstein formula for the diffusion con- 
stant 57, 58]. Different from the total angular momen- 
tum M(r) appearing in Eq. H15|l for the rotating sys- 
tem, Helfand's moment Q{T) of viscosity appearing in 
Eq. H39I) is not a conserved quantity, and this requires a 
different treatment of the canonical distribution for shear 
flow, as will be discussed in Sec. IIV Bl 



7 



B. Canonical Distribution for Sliear Flows 

In the case of Eq. H3t)|l the distribution function /(F) 
is given by 

/(r) - exp (r) - 7Q(r)] } . (42) 

This is the canonical distribution function for shear flow 
in a non-equihbrium steady state. This distribution can 
be attributed to the general form Q of the canonical 
distribution in the case of H{r) = ij('"'=)(r), n ^ 1, 

Ai(r) = g(r) and/ii(r) = 7. 

Now we mention some physical meanings for the 
shear flow canonical distribution function H42() . For this 
purpose we convert the distribution function for 1)42(1 
the canonical variable T into the distribution function 
/"(v(™™',q) for the position q and the velocity v(™°^' 
in the moving frame jf(™°") , and obtain 



/"(v(™°'^),q) 



^ 2 



(mov) 



+ f/(q)+?/'(q) 



(43) 



where the function u'(q) is dcflned by 



N 



(44) 



Here u'{q) can be regarded as a potential correspond- 
ing to the inertial force which pushes particles in the 
direction of the large \qjy\ region, namely in the direc- 
tion of the high speed region of the global current V. In 
other words, this potential u'(q) expresses the effect of 
Bernoulli's theorem in the hydrodynamics. As another 
important point, using the distribution /"(v(™°'"\q) we 
obtain 



JV 

E 



— m 
2 



(mov) 



dNT 



(45) 



for any potential J7(q), where d is the spatial dimension of 
the system. Therefore T can be interpreted as the kinetic 
temperature :59f| . It may be noted that the same relation 
with Eq. H45(l is also derived from the distribution g{T) 
defined by Eq. H13() by interpreting the average X oi X 
as the average under the distribution ^(r). 

The canonical distribution H42I) for shear flow corre- 
sponds to the canonical distribution (|17|l for the rotat- 
ing system, but there is a signiflcant difference between 
the two. Owing to Eq. (g^, the distribution ll^ 
is not stationary in the inertial frame jFl™'^)^ namely 



j2(*"'=)/(r) 7^ 0, whereas the distribution l(T7|l is station- 
ary in time in this frame. (However, note that both the 
distributions H17|) and (|42|) are stationary in time in the 
moving frame namely iL^"""'"^ f{T) = 0.) Phys- 

ically speaking, this difference comes from the fact that 
we need some work to sustain the steady current in the 
shear flow, whereas such work is not necessary in the ro- 
tational system because of the total angular momentum 
conservation law 1)16(1 . However the effect of the work to 
sustain the shear flow is not included in the distribution 
((4211 for the shear flow system. In order to include the 
effect of this work we have to generalize the distribution 
(I42II . and introduce the distribution f{T,t) at time t as 



/(r,t) EE exp{-*L(»^)(t-to)}/(r) 



(46) 



= /(r)exp|-/37Vy dsP^y{T,-s + 2to: 

(47) 

where to is the initial time. Here, to derive Eq. I(47() 
we used the relation gill, iL(*"^)ff(™^) (F) = and 

exp{-iL("-)(t-to)}g(r) = g(r)-v JldsP,y{r,-s + 

2to), and deflned Pxy{T,t) by 

P^y{T, t) EE exp [iL^'^-\t - to)} Pxy{T). (48) 

The distribution 146() corresponds to the distribution (O 
in a general formulation of the non-equilibrium canonical 
distribution approach. It may be noted that the distri- 
bution /(r,t) is normalized, namely J dTf{T,t) — 1, as 
well as J dTf{T) — 1. Using the distribution f{T,t) we 
define the average {X{T))t by 



{x{T))t = J drx{r)f{T,t) 



(49) 



for any function X{T) of F. In the rotating system it is 
easy to check the relation f{T,t) = /(F) and {X{T))t = 
X{T) for any function X{T). 

Noting Eq. ((47(1 . the difference of the distribution 
/(F,t) from the distribution /(F) appears as the factor 
exp{— /37V /^^^ dsPxy{T,—s + 2to)}, and we will discuss 
the relation of this factor to the work needed to sus- 
tain the shear flow in Sec. IIV Dl One may interpret the 
canonical distribution /(F) as a steady distribution func- 
tion in the moving frame j^(™o^)^ b^t in order to know 
about the work to sustain the steady flow we have to in- 
vestigate it from the different frame J^i'^^'^) ^ because the 
work to sustain the steady flow is information given by 
looking at the moving system from the inertial frame. 
Therefore, the canonical distribution /(F) should not be 
regarded as an artiflcial test initial distribution, like in 
other canonical distribution approaches for a linear re- 
sponse theory [s^ l37l| . The information about the work 
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to sustain steady flows is essential to calculate transport 
coefficients such as the viscosity, as will be shown in Sec. 



C. Linear Response Formula for Viscosity 

To calculate the transport coefficient from the non- 
equilibrium canonical distribution approach is beyond 
the purpose of this paper. However many works have 
been devoted for this subject [M El S iH IM IM El , 
so it may be meaningful to mention the consistency of 
the non-equilibrium canonical distribution H47|) with the 
linear response formula for viscosity. 

Using the notation H49|l and the quantity Pxij{T) de- 
fined by Eq. H37|) . the viscosity rj is introduced as 



lim 

7-»0 



1 



(50) 



Using the distribution f{T,t) given by Eq. l|T7|) . the 

viscosity rj is represented as 



77 = /3v / dt (p„j(r,t)p,y(r) 



(eg) 



(51) 



where we introduced the notation {X{T))^'^'^'> as the equi- 
librium average of X(r) for any function X{T), namely 

(x(r))('='?) = s(^9)-i/ drx(r)exp{-/3i?(™^)(r)} with 

the equilibrium partition function 2^"^^^ The derivation 
of Eq. (|51|l is given in Appendix ^ (In the same ap- 
pendixIXlwe also discuss two kinds of nonlinear response 
formulas for {Pxy{T))^ with respect to the shear rate 7, 
one of which can be regarded as a natural generalization 
of Eq. I|51|) .) Here, it is important to note that 



p.y{r) = 0, 



(52) 



at any shear rate 7, as also shown in Appendix^ so that 
we obtain the equation — lim^^o Pxy(r)/'^ = meaning 
that the distribution /(F) does not include information 
about the viscosity. Equation (|51|l is the well known 
linear response formula for viscosity 0|. 

The factor exp{-f3jV J*^dsP^y{r,-s + 2to)} in the 
non-equilibrium canonical distribution H47(l gives the dif- 
ference between the two distributions /(F) and f{T,t), 
and plays an essential role in the derivation of the linear 
response formula (|51|l for viscosity. It may be emphasized 
that this kind of factor can be derived from a different 
approach using the Shod equation 0, E3l- The Shod 
equations expresses the dynamics of the velocity corre- 
sponding to v(™°^^ , and has been used in man y nu merical 
and analytical works on shear flow systems pH E3 • In 
the canonical distribution approach using the Sllod equa- 
tions, the time-evolution of a canonical distribution un- 
der Sllod dynamics is considered, and it leads to the dis- 
tribution evolving a time-integral of the shear stress, like 



the distribution H47|l . In Appendix IBl we discuss briefly 
the relation between the Sllod dynamics approach and 
the Hamiltonian dynamics approach used here. These 
two approaches give the same formula (|51|l for the vis- 
cosity. A difference between this approach and the ap- 
proach discussed in this paper is that the Sllod dynamics 
approach is based on distributions of the type (|13|l . so 
that it does not take into account of the inertial force. 
This make discussions of thermodynamic relations (for 
example, the first law of thermodynamics) rather more 
complicated than the approach used in this paper. It may 
also be noted that the Sllod dynamics is different from 
the dynamics for v(™°") from the Hamiltonian iJ (™°-") (r) 
in the moving frame jF^™""', and in the Sllod dynamics 
approach the distribution corresponding to /(F) is just 
an initial test distribution and cannot be interpreted as a 
steady distribution in the moving frame like the 

Hamiltonian dynamics approach discussed in this paper. 



D. Work Needed to Sustain Shear Flows and the 
House-Keeping Heat 

Now we discuss further the information involved in 
the distribution /(F, t), which the canonical distribution 
/(F) does not have. It is the information about the work 
required to sustain the steady shear flow. 

First, the power {W)t to sustain the shear flow at time 
t is estimated by 



{W)t 



Jt 



di7('"°^)(F,t) 
di 

{P.y), 



(53) 



where i?('"™)(F,i) is defined as #(™°^)(F,t) = 
expiiU'^-^^t - io)}-ff^'"°''Hr)- Here, in order to de- 
rive Eq. lO we used Eqs. (r) = and 

Qffimov) (p^ ^) ^ exp{il{^ne) _ to)}VP^y (F) . 

In steady flow systems, the energy added to the system 
by the work needed to sustain the flow must be elimi- 
nated from the system as heat . This type of heat is 
called the " house-keeping heat" , and its special role has 
been emphasized in the construction of a non-equilibrium 
steady state thermodynamics Fig. QJis a schematic 

illustration of the power {W)t required to sustain the 
shear flow and the house- keeping heat. Now we estimate 
this house-keeping heat from the non-equilibrium canon- 
ical distribution approach. For this purpose we introduce 
the non-equilibrium entropy {S)t as the ensemble aver- 
age of S{T) by the distribution /(F, t). Note that S'(F) is 
defined by Eq. , and was already used as the observ- 
able corresponding to the entropy in Eq. (|23ll . A similar 
kind of entropy to {S)t was used in Refs. jsJ, E3 for a 
different form of the distribution /(F). Using Eq. H47fl 
the entropy {S)t is represented as 
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<W>t 




Td<S>t/dt 



FIG. 1: Schematic illustration of the power {W}t required to 
sustain the shear flow and the house-keeping heat T[d{S}t /dt]. 
In this illustration the power {W)t is represented as the power 
to move the upper boundary of the shearing system. This 
power supplies the energy to sustain the shear flow state, 
which is eliminated from the system as the house-keeping 
heat. 



d{S)t 
dt 



> 



(58) 



This is the expression of the second law of thermody- 
namics in the non-equiUbrium canonical distribution ap- 
proach. The inequahty (|58|) means simply that the shear 
flow system produces a positive house-keeping heat con- 
stantly in time. In this sense, the total entropy produc- 
tion {S)t — {S)to diverges as time t go to infinity, be- 
cause the total amount of heat produced by the steady 
visco-elastic shear flow in the infinite time interval is 
infinite |60| . In other words, the system is kept as a 
non-equilibrium steady state by discharging an amount 
of entropy constantly, which is transferred from the ex- 
ternal work. Therefore the inequality H58|) must be dis- 
tinguished from another type of second law of thermo- 
dynamics meaning that an entropy increases in time and 
approaches to a stable value in a relaxation process. This 
type of the second law of thermodynamics, or the ther- 
modynamical stability condition, will be discussed in Sec. 
IV Bl By combining the inequality H58|l and T > with 
Eqs. lf53)) and (j5H|l we have 



{S)t = - 



I dr/(r)in{/(r,-i + 2to)} 



_ p-t+2to 

= S + PjV dsP,y{T,-s + 2to) 

Jto 

= S-p-fV f ds{P^y)^ (54) 



V 



-7(^.,(r))^ >o. 



(59) 



Namely the averaged power (W)ao needed to sustain the 
shear flow must be positive (or zero). This is one of the 
results, which can be checked by numerical simulation, 
as will actually shown in Sec. IVI Al It may be noted 
that the inequality (|59|l implies the non-negativity of the 
viscosity 77 as a special case. 



Therefore the house-keeping heat T[d{S) 1 1 dt] at time t 
is given by 



THERMODYNAMICS FOR SHEAR FLOW 
SYSTEMS 



- -lV{P.y), (55) 
= {W)t. (56) 

This is the balance equation to show that the power 
needed to sustain the shear flow must be equal to the 
house- keeping heat. Using Eqs. (|5?1|l . (|^ and (|5H|l . 
the house-keeping heat and the ensemble averaged power 
{W)t to sustain the shear flow is connected to the viscos- 
ity Tj as T[d{S)t/dt] = {W)t = V?77' + ^(t')- 
The entropy {S)t satisfies the inequality 

{S)t>{S)t„^S (57) 

at any time t (> <o)- The detail of the derivation of 
the inequality lf77|l is given in Appendix [CI Noting 
that {S)t — {S)to = J*^dt[d{S)t/dt] and assuming that 
d{S)t/dt is time-independent in a steady state, we ob- 
tain 



As discussed in Sec. IIVDI the factor 

exp{7V /j*^ dsPrcy{T,—s + 2to)}, which gives the differ- 
ence between the distribution f{T,t) and the canonical 
distribution /(F), includes the effect of the work needed 
to sustain the shear flow, or the house- keeping heat. 
On the other hand, the expression for the first law 
of thermodynamics proposed by Evans and Hanley 
does not include this effect (Otherwise it must include 
time-dependent terms for the sustaining work and the 
house-keeping heat.). Moreover, Ref. emphasized 
that we must subtract the contribution of the house- 
keeping heat from the entropy in order to obtain an 
expression for the first law of thermodynamics for steady 
states. Equation H54|l implies that the entropy minus 
the contribution from the house-keeping heat is given by 

{s)t^Jl ds[d{s)jds] = {s)t+PiVJl ds{p,y{r)), = s, 

which is the entropy defined through the canonical dis- 
tribution /(r), not through the distribution f{T,t). For 
these reasons, (although it may not be impossible to con- 
struct a non-equilibrium steady state thermodynamics 
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explicitly including the effect of the house-keeping heat) 
in this section we construct a shear flow thermodynamics 
based on the canonical distribution f{T) excluding the 
effect of the house-keeping heat, and show that it is 
consistent with the Evans-Hanley thermodynamics. 



A. First Law of Thermodynamics for Shear Flows 

As the first thermodynamic property of the shear flow 
system, we consider the first law of thermodynamics. Its 
derivation is quite similar to that for the rotating system 



discussed in Sec. IIII Bl so it is given rather briefly. 

In the shear flow system, the entropy S is given by Eq. 
(ESI or 



S = lnE + P H(.ine)_^Q 



(60) 



using the canonical distribution (|42|l . As in the rotat- 
ing system, the free energy p(."^°^) in the moving frame 
j7(mov) jg defined by Eq. and is connected to the 

partition function S by Eq. (|25(l . (Here, it may be 
noted that the free energy 

pimov) 

can also be expressed 
as pi^ov) ^ (i7(™o'')(r))t - T{S)t using the averaged 
energy {H^"'°'"'> {T))t and entropy {S)t related to the dis- 
tribution /(r,t) which includes information about the 
house- keeping heat.) Using this free energy we 
obtain Eqs. and 



ainS 



-f3Q 



(61) 



which are equivalent to dlPF^"^ 
Therefore we obtain 



(62) 



noting Eq. Using the relations F'-'""'^ = p(niov) _ 

pimov) _,_ pg g^^^ J^U^ ^ p(^ne) ^ p-g 



7(5, TJC™"") 

also obtain 



dF^'"""^ = -SdT + jdQ, 



dHi^"-") = TdS - Qdj, 



dH(^ne) ^ pdg _,_ 



(63) 



(64) 



(65) 



The two energies 
frames are connected 



and iJl™"") in the different 
by a Legendre transformation, 
namely = i7(™o") _ ^ dlK"^"") /dj\-g, as well as 

the two free energies in the different frames. It is clear 
that the Evans-Hanley expression ^ for the first law 
of shear fiow thermodynamics is the relation H64|l where 
£ ^ and f = -Q. 



B. Thermodynamic Stability Conditions for Shear 
Flows 



As the second thermodynamical property of the shear 
flow system, although we omitted to discuss it in Sec. 
IIII Bl for the rotating system, we consider a stability con- 
dition for shear flow f6l| . 

We consider a small part A of the macroscopic shear 
flow system, in which averages of energy, entropy and 
Helfand's moment of viscosity in the inertial frame J^i™'^) 
are given by S and Q, respectively. The other 

part TZ of the system, which is much bigger than the 
system A and is called the "environment" or "reser- 
voir", has the thermodynamical values Tq, 5*0 and 70 
of the temperature, the entropy and the shear rate, re- 
spectively. Now, we consider moving an inflnitesimal 
amount of energy as heat —TodSo from the reservoir 
TZ into the small system A. In this process the total 
entropy must increase: dS + dSo > 0. By combin- 
ing this inequality with the flrst law of thermodynamics 

H65|) . we have 
7oQ) 



dH(^ 



- -TodSo 

rf7j(»ne) _ p^ds _ j^dQ = d{H("^<=^ - TqS - 7oQ) < 0, 
using the fact that the reservoir TZ is so big that Tg and 
7o do not change in this process. This inequality means 
that the quantity H'^"^'^'> — TqS — joQ always decreases 
and reaches a minimum at a stable state. In other words, 
if we force a change to the values of (*"'^) , S and Q at 
the stable point by SH^^^'^^ , 5S and 6Q, respectively, then 
the inequality SH('"'') - TqSS - jqSQ > must be satis- 
fled as the stability condition for the shear flow system. 
This simply leads to 



^odQ based on Eq. 

-ToS 



> 



(66) 



for any inflnitesimal deviations SS and 5Q. By a well 
known technique used in thermodynamics (See, for ex- 
ample, Ref. or Appendix IdI). the condition H66|l is 
equivalent to 



dS_ 
df 

dQ 



>0, 



> 0. 



(67) 



(68) 



The condition (j67|l simply means that the speciflc heat at 
constant Q is always positive at a positive temperature 
T. To understand the condition H68|l we note 



dQ_ 



(69) 



as shown in Appendix^ Therefore, combining Eq. H69|l 
with the inequality H68() we obtain 
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> 0. 



(70) 



Namely, the stability condition (|69|) means the positivity 
of the correlation function for Helfand's moment Q(r) of 
viscosity. 

Based on Eq. Evans and Hanley claimed the in- 
equality 9£/97|t > as a stability condition for shear 
flows ^s, zy| . This inequality is incompatible with the 
inequality in the case of ^ — —Q. This difference 
comes basically from the fact that they discussed a ther- 
modynamic stability condition using the energy H^"^°^\ 
whereas we discussed it using the energy Obvi- 
ously the correlation function of Q cannot be negative be- 
cause of — Q = {Q — QY ^ 0: so noting Eq. (jSHJ we 
cannot justify the stability condition claimed by Evans 
and Hanley in the canonical distribution approach. 



C. Relations Between Canonical Averages 

So far we have introduced two types of canonical aver- 
age X and {X)t, and in Sec. IVII we introduce the usual 
time average. It is very important to distinguish between 
these averages. The thermodynamic relations discussed 
in Sees. IV Al and IV Bl are the relations for the ensemble 
average X of observable X{T) using the canonical dis- 
tribution f{T). On the other hand in numerical simula- 
tions using the Shod equations with an isokinetic thermo- 
stat (as in Sec. IVI|I . the values obtained are the mixed 
ensemble-time average {X)t for the distribution f{T,t) 
in the limit i — > oo. Therefore it is important to obtain 
an explicit relation between these two different ensemble 
averages. 

For any function X{r) the relation between the two 
ensemble averages X and {X)oo is 



= X - /37V 



dt 



X{T,t)~X P^y{T)-P, 



(71) 



where X{T,t) = exp{iL(^"^)(t - to)}X{T). The deriva- 
tion of Eq. H71|l is given in Appendix IeI A similar equa- 
tion for the canonical distribution approach using the 
Sllod equations is shown in Ref. [sof . 

From relation ifTTl) . if the fluctuation X{T,t) — X of 
X is weakly correlated to the shear stress Pxy{T), then 
the ensemble average X can be nicely approximated by 
the average {X)oo- However one must notice that the 
justification for such an approximation strongly depends 
on the quantity X we consider. A typical example is the 
case oi X — Pxy{T), in which we must not neglect the 
second term on the right-hand side of Eq. H71|l . because 
in this case the first term on the right-hand side of Eq. 
(|71|l is zero, namely Pxy — 0, as shown in Appendix 1X1 



One should also notice that the second term in Eq. H71|) 
is small near equilibrium, because it includes the non- 
equilibrium parameter 7 as a factor. 



VI. NUMERICAL SIMULATIONS OF SHEAR 
FLOW 

In this section we show numerical results for some 
quantities which have appeared in the preceding sections 
II VI and Ivl and check the results obtained there. 

For this numerical calculation we use a two- 
dimensional square system consisting of N particles with 
a square shape and side length L{= VV). The particle- 
particle interaction is given by the isotropic soft-core 
pair-potential 




) m r < ro 
in r > vq 



(72) 



with a positive constants vq = 1.5 and k = 1. The par- 
ticle number density p = N/V is 0.8. The mass m of 
the particle and the kinetic temperature T are both cho- 
sen as 1. The number of particles is iV = 450, except 
in Sec. IVIDI in which the iV-dependence of a quantity 
will be discussed. We use the Sllod equations with Lees- 
Edwards boundary conditions and the isokinetic thermo- 
stat so that the kinetic temperature (given by Eq. (|45|l 'l 
is kept constant 17]. (This dynamics is explained in Ap- 
pendix ^ more explicitly.) A predictor-corrector method 
62] of 4-th order is used to carry out these numerical sim- 
ulations with time step of At = 0.001. In this algorithm 
the sum of the "thermal momentum" = pj — rnV{q^j) 
is zero in the both coordinate directions. 

We use the notation {X) for the time-averaged value of 
any quantity X given by this numerical simulation. To 
calculate this average we used data over more than 10^ 
time steps omitting the first 10'* time steps. (We checked 
that 10"* time steps is much longer than the relaxation 
time of the time-correlation function for the thermal mo- 
mentum.) This should correspond to the ensemble av- 
erage {X)oo used so far. This is supposed by the fact 
that we can calculate the viscosity from the time-average 
(Pxy) in this simulation, based on Eq. (|50|l assuming 
{Pxy) 00 = {Pxy)- We calculate 7-dependences of three 
quantities: {Pxy), {Q) and (Q^). We use {Pxy) to discuss 
the power to sustain the flow and the house-keeping heat 
given by Eq. ifSsl) . The quantities {Q) and (Q^) are used 
to discuss the behavior of Helfand's moment of viscosity 
and the thermodynamic stability condition Ht)8|) . Here, it 
is assumed that the quantities Q and Pxy are not strongly 
correlated with each other because of the relation H41|l . 
and in the case oi X = Q ot the second term of the 
right-hand side of Eq. (|71() may be small compared to 
the first term in the small shear rate case. This implies 
that the behavior of the time-averages of Helfand's mo- 
ment of viscosity and its correlation function are not so 

different from the ones for Q and — Q , respectively, 
near equilibrium. 
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A. Work Needed to Sustain the Shear Flow 

The first numerical result is for the power W /V per 
unit volume (V = 562.499 • • • ) for the work required to 
sustain the shear flow, or equivalently in a quantitative 
sense, the house-keeping heat_ per unit volume. It is 
given as the time-average of W /V = —jPxy{T), based 
on Eqs. H55|l and (|56|l . Figure |21 shows the shear rate 7 
dependence of the time-averaged power {W)/V per vol- 
ume. The inset is the same graph except showing it in a 
wider shear rate region. (Note that we used a linear-log 
scale in this inset, whereas we use a linear- linear scale 
for the main figure.) Following from the inequality (|59|l . 
the power to sustain the flow shown in Fig. |21is always 
positive (or zero). 

y = apX2 y = cip'x 




1 2 3 4 5 6 



Y 

FIG. 2: The averaged power (W) /V divided by volume V 
needed to sustain the shear flow, as a function of shear rate 7 
as a linear-linear plot. The solid line is the fit to a quadratic 
function and the dashed fine is the fit to a linear function 
(vafid for 7 > 3). Inset: the same graph as a linear-log plot 
including a wider range of 7. 

It may be noted that the averaged power (W) needed 
to sustain the shear flow should be an even function of 
7, because it should be invariant under a change of sign 
of the shear rate 7. In Fig. |21 we fitted the numerical 
data to a quadratic function y = apx^ with the fitting 
parameter ap = 1.51309. Near equilibrium 7 < 0.5, the 
graph is nicely fitted by this quadratic function. 

As the shear rate increases, a region in which the value 
of {Pxy{T)) is almost independent of 7, (namely the re- 
gion fitted by a linear function y = with the fitting 
parameter a'p = 2.00005), appears Mi, and after that 
the region the string phase appears 0, ^| (the re- 
gion 7 > 8 approximately in the inset to Fig. The 
string phase can be checked not only by the string-type 
arrangement of particle positions but also by strong time- 
oscillating behavior of the time-correlation functions for 
quantities such as the potential energy, the shear stress 
and so on fTS|. 

For isokinetic thermostat, the house- keeping heat is 
also given as the time-average of the thermostat term 



(the term —a{t)p'j{t) explained in Appendix IB)i . We 
checked numerically that this quantity is equal to the 
time-average of —jPxy{T). 

B. Helfand's Moment of Viscosity 

As a second numerical result. Fig. 01 shows the graph 
of the time-average (Q)/V of the Helfand's moment of 
viscosity per unit volume V as a function of shear rate 7. 
It (almost) takes the value at the equilibrium 7 = 0, 
and increases linearly as a function of 7. In this figure 
we also give a fit to a linear function y = Q!q7, with the 
parameter value aq = 150.043. 



y = aox 
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FIG. 3: Time-average (Q) /V of Helfand's moment of viscosity 
divided by volume V as a function of shear rate 7. The solid 
line is the fit to a linear function. Inset: the same graph 
including a wider range of 7. 

To explain this linear behavior for Helfand's moment 
of viscosity as a function of shear rate, we simply note 
that 

N N 

{Q)^J2(^^vPj-^+^J2(i%) (73) 

with the x-component pjx = Pjx — IQjy of the ther- 
mal momentum of the j-th particle. Our numerical 
calculations show that the value of the first term of 
the right-hand side of Eq. H73fl is extremely small (or 
zero) compared to the value of its second term, namely 
'Y^^=i{ljyPjx) ~ 0. Moreover, using a homogeneous 
continuum assumption for the fluid, the value of the 
quantity '^jLii^jy) appearing in the second term of 
the right-hand side of Eq. H73|) can be estimated as 
Y.^=i{<fjy) ~ NL-^J^dyy'^ = These estima- 

tions lead to {Q)/V ~ {N/3)j = I5G7, which explains 
the value of the fitting parameter ag. 

The time-averaged Helfand's moment (Q) of viscosity 
should be at least an odd function of shear rate 7, be- 
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{Q') {Qf 



cause the infinitesimal deviation ^dQ giving the energy region. This point may be explained by noting 
change dH^™'^^ in the inertial frame by Eq. Ht)5|) must be 
invariant under the change of sign of the shear rate. It 
may be noted that this linear dependence for the time- 
average of Helfand's moment Q of viscosity with respect 
to shear rate is satisfied not only in the near equilibrium 
region but also even in the string phase region, shown in 
the inset to Fig. O possibly because the Shod equations 
are a homogenous shear algorithm. 

It may be noted that in our simulations Helfand's mo- 
ment of viscosity can be generally changed discontinu- 
ously in time, when a particle steps over a boundary in 
the direction orthogonal to the global shear flow. How- 
ever it should be a small boundary effect which can 
be neglected in the thermodynamic limit N ~* oo and 
p = const, and our numerical calculations gave a good 
convergence for the long time-average of Helfand's mo- 
ment of viscosity. 



N N 

X! X! iljyIkyPjxPkx) 

]=1 k=l 

N N 

j=i k=i 

(74) 

using the thermal momentum component pj^- Here 
we assumed that the time-average of the linearly de- 
pendent terms for the thermal momentum can be ne- 
glected. As the two terms I]j=iEfc=i( 

QjyQkyPjxPkx/ 

^ ^ '"^ ^ can be considered 7- 
2\ - (Q)2 can be fitted 



C. Correlation Function for Helfand's Moment of 
Viscosity 

As the last example, Fig. ^ shows the shear rate de- 
pendence of the correlation function ((Q^) — {Q)'^)/V of 
Helfand's moment of viscosity divided by the volume V. 
This figure shows that this correlation function is always 
positive at least for 7 < 10, following the thermodynamic 
stability condition (jSHl- 




FIG. 4: Correlation function ((Q^) - {Qf)/V of Helfand's 
moment of viscosity divided by volume V as a function of 
shear rate 7 in a linear- linear plot. The solid line is the fit to 
a quadratic function. The inset: the same graph except for 
that it includes a wider region of 7 and is a linear- log plot. 



The correlation function for Helfand's moment of vis- 
cosity should be an even function of the shear rate. 
Noting this point, in the small shear rate region of 
Fig. 2] we give the fit of numerical data to the func- 
tion y = Q!(7o + ciC2X^ with the fitting parameter values 
otco = 46.749 and ac2 = 1349.38. The graph is nicely 
fitted by this quadratic function in the small shear rate 



and Ei=i Efe=i Ej=i {TjyTky 
independent, the correlation {Q 
by a quadratic function of 7. 

In the inset to Fig. 2) we give the shear rate depen- 
dence of the time-average of the correlation function for 
Helfand's moment of viscosity per unit volume V in a 
much wider region of shear rate in a linear-log scale. 
(Note that we used a linear-linear scale in the main fig- 
ure of Fig. 21) It should be noted that a rapid drop of 
the value of this correlation function occurs in the string 
phase region. In the intermediate region, which is ap- 
proximately the region 2.5 < 7 < 8 in Fig. 01 between 
the region fitted by the quadratic function of 7 and the 
string phase region, fluctuations in the value (Q^) — {Q)'^ 
become much larger than in the other regions, and their 
values in Fig. ^are less reliable. 



D. Remarks in Connection with the Isokinetic 
Thermostat Dynamics and the Canonical 
Distribution Approach 

Shod dynamics with the isokinetic thermostat used 
in this section has been used very frequently to simu- 
late shear flows. It is supposed to reproduce the value 
of shear stress predicted by a canonical distribution ap- 
proach and succeeded even to reproduce some real 
experimental values '12^1. However one must notice that 
strictly speaking the time-average from Sllod dynamics 
with the isokinetic thermostat does not always reproduce 
the ensemble average for the non-equilibrium canonical 
distribution used in this paper, even in the equilibrium 
state where 7 = after taking the thermodynamic limit 
N 00 (and p — const). Now we discuss a couple of 
examples illustrating these ensemble differences. 

First, in the numerical simulations used in this section 
the sum of the thermal momentum = (pjxiPjy) = 
Pj — mV(qj) over the particle number j in each direction 
is zero at all time, meaning that there is a constraint on 
the values of the thermal momenta, that is X^jLi-Pja; ^ 0- 
On the other hand, in the canonical distribution ap- 
proach all components of momenta can be treated as in- 
dependent variables. This difference, for example, causes 



the different averaged values for J27=iJ2k=iiPjxPk 
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and J2j'=iJ2k=iPr-cPkx- Actually the value of 



'Ej=iT,k=i(PjxPkx) = {iT,j=iPjx){T,k=iPkx)) is zero 
as each bracketed sum is individually zero. The value 



of J2k=iPjxPkx however, is given by mNT in the 

canonical distribution because of pjxPkx = rnTSjk- 

Second, the isokinetic thermostat used in the simula- 
tions of this section keeps the kinetic energy constant 
so that the distribution function for the kinetic energy 
is given by a delta function. This is different from the 
distribution of kinetic energy derived from the canonical 
distribution, because there is always a non-zero fluctua- 
tion of the kinetic energy around its mean value in the 
canonical distribution. Refs. (63.. .64] tried to modify the 
canonical distribution to give consistency with the isoki- 
netic thermostat, but it is not obvious that we can justify 
the shear flow thermodynamics based on such a modified 
canonical distribution. 




FIG. 5: Time-averaged quantity $i 



a function of particle number A*' for Sllod dynamics with an 
isokinetic thermostat in a square at equilibrium 7 = with 
the particle density p — 0.8 as a log-linear plot. Here "i/xy and 
are defined by '^^y = J^f^i J2k=i{ljy<lkyPjxPkx)\-,=o and 
^xy = "iT ^^^-^ 1-^=0, respectively. The length of error 
bars in this figure is given by 2|"I>a;y — ^yx\- The solid line is 
the value 3/4, which is explained in the text. 



As a concrete example of these ensemble dif- 
ferences, let's consider the first term '^^y = 
Y^f=iY^k=ii9jv1kyPjxPkx)U=Q appearing on the right- 
hand side of Eq. H74() at equilibrium 7 = 0. Assuming 
that the variables qjy and pjx are independent, and that 
{ijyQky) only depends upon whether j = k or j k, then 

N N 

i<ljy1kyPjxPkx) 

3 = 1 k=l 

N N N 

~ (liy) Y (P%) + iliyl^y) Y ^JxPkx) (75) 

j=i j=i fc=i(fc#j) 

If qiy is uniformly distributed between and L then 
{qly) « and {qiyq2y) ~ L^/4. The difference be- 



tween 'i!xy and 'i!xy can be calculated as pjxPkx — 
in j ^ k for the canonical average, whereas for the 

time average Y^j^kiPjxPkx) = -J2f=i{P%)- I* follows 

that ~ mNTL'^/3 and ^^y « mNTL'^/12, so 



-xy 

EE \^xy - '^xyl/'^xy ~ 3/4. Figure is the graph 
of the normalized difference ^xy as a function of system 
size N for square systems at fixed density p — 0.8 in 
numerical simulations. The length of error bars in this 
figure is given by 2|$a;j, — ^yx\ which must be zero in the 
square cases. Figure suggests that ^xy is in excellent 
agreement with the value of 3/4 given above. 



VII. CONCLUSION AND REMARKS 

In this paper we have discussed a canonical distribu- 
tion approach to non-equilibrium steady flows for the 
purpose of constructing a steady state thermodynam- 
ics from solid statistical mechanical foundations. Using 
the Lagrangian technique of classical mechanics we in- 
troduced the energy in the moving frame by separating 
the velocity of the global steady flow. A canonical distri- 
bution based on this internal energy was introduced. As 
one application of this distribution, we showed that the 
well known thermodynamics of rotating systems can be 
derived from this canonical distribution. Our special con- 
cern in this canonical distribution approach was steady 
shear flows and their thermodynamics. Evans and Han- 
ley proposed a first law of thermodynamics of the form 
d£ = TdS — Qd'y relating energy £, temperature T, en- 
tropy S and shear rate 7. Here we derived this shear 
flow thermodynamics based on our canonical distribu- 
tion approach, and showed that the quantity Q is given 
by the average of Helfand's moment of viscosity, the tem- 
perature T is the kinetic temperature derived from the 
thermal kinetic energy, and £ can be interpreted as an 
internal energy. The roles of the work required to sus- 
tain the shear flow and the heat removed to compensate 
it (the house-keeping heat) was emphasized in the justifi- 
cation of the linear response formula for viscosity, which 
is derived from our shear flow canonical distribution ap- 
proach. We introduced a non-equilibrium entropy, and 
showed that it increases in time and the house-keeping 
heat based on this entropy has the same magnitude as 
the power needed to sustain the steady flow. This dis- 
cussion led to the non-negativity of average of —^Pxy 
with the shear stress Pxyi meaning that the power needed 
to sustain the shear flow and the house-keeping heat is 
always non-negative. We discussed the thermodynamic 
stability condition for the shear flows, one of which is 
equivalent to the positivity of the correlation function 
of Helfand's moment of viscosity. Our results were in- 
vestigated in numerical simulations of two-dimensional 
many-particle systems with soft-core interactions, whose 
dynamics is determined by the Sllod equations with an 
isokinetic thermostat. 

To construct the canonical distribution for shear flow, 
we used the analogy of shear flows and rotating systems. 
These two systems are steady flows whose magnitude is 



15 



proportional to a component of position vector: the dis- 
tance from the rotating axis in the rotating system, and 
the position component orthogonal to the flow in the 
shear system. Both systems also have clear parameters to 
characterize their currents: the angular velocity in the ro- 
tating system and the shear rate in the shear flow. On the 
other hand, we also emphasized some differences between 
these two systems. The biggest difference may be that 
the total angular momentum in the canonical distribu- 
tion of the rotating system is time-independent, whereas 
Helfand's moment of viscosity appearing the canonical 
distribution of the shear flow is not constant. This led 
to the necessity to consider the work needed to sustain 
the steady flow and the house-keeping heat in the shear 
flow system, and plays an essential role in the derivation 
of the response formula for viscosity. 

One may easily notice that the canonical distribution 
approach discussed in this paper can be generalized to 
more general steady flows than the rotating system and 
the shear flow system. One of the restrictions in our 
canonical distribution approach is that we have to know 
the global velocity distribution V a priori. In this sense 
this approach is not appropriate to determine the global 
velocity distribution under some external constraints, 
etc. It is also crucial that we know a priori an external 
parameter that specifies the amount of the global flow, 
like the angular velocity or the the shear rate. This pa- 
rameter is treated as a thermodynamic quantity in the 
expression for the first law of thermodynamics. 



<Pxx(r')> o 



<Pyy(r')> A 




FIG. 6: Time-averages of the diagonal components {Pxx(r')) 
(circles) and {Pyy(T')) (triangles) of the pressure tensor with 
the thermal phase space vector T' as functions of shear rate 
7 in SUod dynamics with an isokinetic thermostat in a linear- 
linear plot. The inset: the same graph except for including a 
wider region of shear rate and in a linear- log plot. 



An important future problem in shear flow thermo- 
dynamics based on the canonical distribution approach 
is to discuss the pressure in this framework. Refs. 
[27I [2s. 30] introduced the pressure V simply by 
adding the term —VdV on the right-hand side of Eq. 



(Q. For this term it was conjectured that the pressure V 
would be equal to the minimum eigenvalue of the pres- 
sure tensor [s^l- However one should notice that non- 
equilibrium systems such as the shear flow system are 
not generally isotropic, so that the pressure defined by 
—d£/dV may depend on which direction we change the 
volume V. Actually, as shown in Fig. El the numerical 
simulations using Shod dynamics in Sec. IVII show that 
the time averages of PxxiX') ^'^'^ Pyyi^') different 
from each other at non-zero shear rate. (Here F' is the 
"thermal phase space vector" introduced as the vector 
in which the momentum pj in the phase space vector is 
replaced by the thermal momentum pj = pj — mV(qj).) 
Noting that usually the pressure is calculated by the 
arithmetic average of these time-av erag es (or ensemble 
averages), (See Ref. (T^, also Ref. "Sffl for its justifica- 
tion using the microcanonical distribution.), this suggests 
that if the pressures in the x and the y-directions are 
given by averages of Pxx(r') and Pyy{T'), respectively, 
then the pressure is direction-dependent in shear flow sys- 
tems. The quantity {Pxx(r')) — {Pyy(^')) is called the 
" normal stress" and a non-zero value is one of the impor- 
tant properties of visco-elastic fluids 0, 0, |6l| . There- 
fore it is important to understand whether such a prop- 
erty is compatible with the thermodynamical framework 
discussed in this paper, in other words, to discuss the first 
law of thermodynamics in which the averages PxxiX') 
and Pyy{T') are included as the x and j/-components of 
the pressure, respectively. It may be noted that a simi- 
lar question can be asked for rotating systems. We leave 
discussion of these points for the future. 

As mentioned in Sec. IV CI the thermodynamic rela- 
tions, (|65|l and (|70|) derived in this paper, are relations for 
the ensemble average (|19f) under the canonical distribu- 
tion /(r). On the other hand the numerical calculations 
discussed in Sec. IVII give the average (|49|l under the dis- 
tribution /(r,oo). Although these two averages are re- 
lated by Eq. H71|l . it is still an open question to calculate 
the canonical average (|19|l . required for the thermody- 
namic relations, from the dynamical evolved canonical 
average (I49|l in numerical calculations. 

Originally, Evans and Hanley introduced their shear 
flow thermodynamics to discuss non-analytical properties 
of the pressure, viscosity and the internal energy as func- 
tions of the shear rate. Such non-analytical properties 
are predicted by mode-coupling theory 38, 65|, and are 
supported by some numerical calculations ,11, .66. 
However, recently some numerical works suggest that 
the shear rate dependence of the pressure is rather ana- 
lytic near equilibrium, except at the triple point 12, 13| . 
Moreover, even at the triple point the non-analytic de- 
pendence of the pressure is not completely convincing 
M. It may also be noted that some theories, that pre- 
dict an analytic dependence of the pressure and the vis- 
cosity with respect to the shear rate, have been proposed 
0, Hoi [63 ■ In this sense it is still an interesting problem 
to discuss shear rate dependences of the pressure, the 
viscosity and so on using shear flow thermodynamics. 

There are also questions about the numerical Simula- 
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FIG. 7: a:-component {Tx{T')) (circles) and y-component 
{Ty(T')) (triangles) of the Icinetic temperature with the ther- 
mal phase space vector T' as functions of shear rate 7 in the 
Sllod dynamics with the isokinetic thermostat. The inset: the 
same graph except for including a wider range of shear rate. 



tions of shear flow themselves, from the point of view of 
the canonical distribution approach. Some such prob- 
lems were already mentioned in Sec. IVI Dl As an- 
other potential problem, we mention the direction de- 
pendence of the thermal kinetic energy. To discuss this 
point we introduce the quantities Tx(T') and Ty{T') as 

rfe(r') = (2/A^)X;j^ip|fc/(2TO) with the thermal mo- 
mentum component pjk = Pjk—rnVkicij), where 14 is the 
fc-component of the global current density V. The arith- 
metic average of ensemble averages of {Tk{T')}k over the 
component k gives the kinetic temperature, so we may 
interpret the quantity Tk{T') as the observable for the 
" /c-component of the temperature" . The canonical dis- 
tribution approach discussed in this paper claims that 

the ensemble average Tk{T') of the quantity rfc(r') is k- 
independent, in other words the kinetic temperature is 
direction- independent, although we should note a differ- 
ence in the two averages Tk{T') and {Tk{T'))oc- Figure 
13 shows the graphs of (T^(r')) and (Ty{T')) as func- 
tions of shear rate 7 from numerical simulations using 
the Sllod dynamics with an isokinetic thermostat, used in 
Sec. I VII This figure shows that the kinetic temperature is 
direction-dependent at least in large shear rate cases. As 
a topic related to this point, it may be noted that in the 
isokinetic thermostat the heat is removed from any com- 
ponent of kinetic energy of any particle uniformly. This 
gives great simplification in the formula and numerical 
calculations and keeps a similar dynamical structure to 
Hamiltonian dynamics leading to the numerical observa- 
tion of the conjugate pairing rule for the Lyapunov spec- 
trum [69l |. but its physical justification as a mechanical 
thermostat is not completely convincing. For example, 
one may use other types of thermostats in which the heat 
is removed from the particles near the walls or from the 
kinetic energy component orthogonal to the walls (70Ll7ll| . 



These different thermostats might give, for example, dif- 
ferent values of {Tx{T')) and (Ty(r')) from the isokinetic 
thermostat. To check the shear flow thermodynamics for 
such types of thermostat remains an open problem. 
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APPENDIX A: RESPONSE FORMULA FOR THE 
VISCOSITY FROM THE CANONICAL 
DISTRIBUTION APPROACH 



In this appendix we give a derivation of the linear re- 
sponse formula H51|l for viscosity from the definition l|5U|l . 
as well as a derivation of Eq. H52|l . We also discuss the 
two kinds of nonlinear response formulas for {Pxy{T)) ^ 
with respect to the shear rate 7, one of which is a simple 
generalization of the formula 1)51(1 . 

First we note that the partition function S can be 
rewritten as 



S = y" drexp. 


'-/377(™°'')(r)} 


= J drexp- 


' -iL^^'^^t - to) 


X exp 1 


_^^('no«)(r)| 


= J drexp- 





X exp <j -/37V / dsP^y{r, -s + 2to) } (Al) 

J to 



where we used the relations exp{—iL^™'^\t — to)}l — 1, 



^j^{me) jneaning to take its Hcrmitian 



conjugate) and a similar derivation to that in Eq. 147|l . 
Equation l(Al|l means that both the distributions /(F) 
and /(r, t) are normalized with the same partition func- 
tion S. The partition function S given by Eq. (|A1() must 
be time- independent, so that we obtain 
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= -:- 



as 

dt 



/37V / dsP,y{T,-s + 2to] 

to 



= -f3jV J drP,y{T,-t + 2to) 
xexp{-/3i7(™°'')(r)} 
X exp 

= -/37VS J dTP,y{r,~t + 2to) 

xexp{-*L(™^)(i-io)}/(r) 
= -/37VS J drexp|-iL(™")(i-to)} 

xp, ,(r)/(r ) 



(A2) 



noting the definition of tlie average (dnj. From Eq. (|X2ll 
we obtain Eq. (|52|1 . impiying that the viscosity calculated 
from the canonical distribution /(F) is zero. 



On the other hand, using the average (|49|) by the dis- 
tribution f{T,t) given by Eq. H47I) we have 



(P,y(r))^ = lim / drP,j,(r)/(r)exp -/37V / dsP^y{T,-s + 2to) 



dTP^y{r)f{r) - /37V J dsj drp^y(T)f{T)p,y{r, -s + 2to) + o{Y) 

P^-PlV ds J dTP,y{r)f^''^\r)P^y{T,-s + 2to)+0{^^) 

-hV ds I dTP,y{T)f^''^\r)exp{-tL^'"-\s-to)}p„j{r)+0{^^) 

-/37V J ds j dT [exp{iL(™^)(s-to)}P.y(r)/('=')(r)] P.,y{T) ^ Oi^"") 
r°° / ~ \ (eg) 

-PjV dslP,y{r,s)P,yiT)\ +0(7') 

J to 



r 



(A3) 



with the notation /'^'^^^(r) = lim-y^o/(r), where we 
used Eq. the relation (ii(™'=))t -iif™'^), and 

exp{iL^""'Hs-to)}f^'"i\r) = f'^^'^HT). Equation llall 
leads to the linear response formula (|51|1 for viscosity. 
Next, using Eq. 147|l we have 



The solution of the time-differential equation (jA4|) of the 
function /(F, t) with the initial condition 



/(r,to) = /(r) 



(A5) 



df{T,t) 
dt 



f3"fVP.^y{r, -t + 2io)/(r, t) (A4) is represented as 



/(r, t) = /(r) + (-/37V)" / dsi r ds2 r ds^-- - r ' dsn 

xP,y{T, -Si + 2to)P^y{T, -S2 + 2to) ■ ■ ■ P^y{T, -.s„ + 2<o)/(r). 



(A6) 



From Eqs. H52() and (|A6|) we derive 
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C>0 nOO pSi pS2 rs^ 



ds„ 



xP,y{T)P,y{r, -Si + 2to)P^y{T, ~S2 + 2to) ■ ■ ■ P:,y{r, - 



2to). 



(A7) 



This expresses a nonlinear response formula for an aver- 
age of the shear stress Pxy{T) with respect to the shear 
rate 7 in the form of its multiple time-correlation func- 
tion. The formula H51|) can be derived directly from Eq. 
(pTTj) . using the relations X{T)\^=q = {X{r)\^=o)'-'"''^ in 
any function X{T) of T and (iL(™^))t = It 
may be noted that the multi-time integral functions on 
the right-hand side of Eq. (|A7I) can be 7-dependent 
because of the 7-dependence of the function /(F), so 
strictly speaking Eq. HA7|) is not an expansion formula 
for {Pxy{T))^ with respect to the shear rate 7. 

It may be meaningful to show another type of nonlinear 
response formula for the average of the quantity P^yiT) 
with respect to the shear rate 7, using a Green's function 
G defined by 



G= lim 



(AS) 



For this purpose, first we note a formal identity 



lim exp(iL("^)A 
= lim^e j dtexp{-ei} • exp|iX(*"^)i| 



= lim e 



and e[e 

. iL(™'=)] = 1 - 



(A9) 



Equation 



= i-gl(™"), 

where we used Eqs. 

(|A9(I is an analogous technique to that used in quantum 
scattering theory j72j in which the Hamiltonian operator 
instead of the Liouville operator L(™^) is used. Using 

Eq. (jsnj we have 

= lim /rfr/(r)exp|iL(™^)(t-to)|P.j;(r) 



dTf{T) 



1 - gl(™'=) 



n=0 



n 



1 - GL^"'"^ 



(eg) 



(AlO) 

This is the formula which we wanted to derive. It may be 



noted that the quantity ([Q(r)]"[l - G'i(™'=)]F^y(r))('=«) 



appearing on the right-hand side of Eq. (|A10p is 7- 
independent, so Eq. lAlOp can be regarded as a real 
expansion of {Pxy{T))^ with respect to the shear rate 
7, different from the formula (|A7|) . Another merit of the 
formula (|A10|) is that wc do not have to calculate a time- 
integral in the interval [0,oo], which is required in the 
formula ljA7p . As a special case of the formula ljA10|l . 
using Eq. (|50ll and the fact that zero-th order of the 
quantity (^^^^(r))^ must be zero, we obtain 



1 - GL^'""^ 



Pxyi^, 



(eq) 



0, 



/3(Q(r) 



1 - Gi 



(eq) 



(All) 



(A12) 



Eq. (|A12|I is another type of the linear response formula 
for the viscosity. 



APPENDIX B: CANONICAL DISTRIBUTION 
APPROACH USING THE SLLOD EQUATION 

In this appendix we discuss briefly the non-equilibrium 
canonical distribution approach using the Shod equation 
for shear flows, and the relation of this approach to that 
discussed in the text of this paper. 

The Sllod equation of shear flows is introduced as the 
equations for d-dimensional vectors <ij{t) and Pj{t) ^3 



dt 



m ■' 



dv',(i) _ 9t/(q^(t)) 
dt dq'At) 



(Bl) 



(B2) 



where C/(q'(t)) is the potential energy as a function of the 
position q'(t) = {^[{t) , q'2{t) , ■ ■ ■ ,q^(t)), and the d x d 
matrix Q = (Qjk) is defined by 



e 



jk 



1 in J = 1 and k 
otherwise 



(B3) 



for j ~ 1,2, ■ ■ ■ ,d and fc = 1, 2, • • • , d. 

The second term on the right-hand side of Eq. IjBll) is 
added so that the quantity p' = {p[, P2, • • • , p^) means 
the mass times the velocity in the moving frame by the 
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velocity of the global shear flow, the so called thermal 
momentum. The second term on the right-hand side 
of Eq. ljB2|) is added so that we derive the equation 



md^q'Jt)/de ^ ^dU{q'{t))/dq'Jt) of the position q'(t) 



from Eqs. and 6^ = 0. Therefore the pa- 

rameter 7 dependence does not appear in the equation 
of q'(t) only. 

Usually the Sllod equation is used with a thermostat, 
such as the isokinetic thermostat (17i |. in order to make 
a model for the system driven by a shear rate with an 
attached heat reservoir which removes (as the house- 
keeping heat) the energy generated inside the system by 
the shear and maintains the temperature of the system 
constant in time. The isokinetic thermostat is expressed 
by the term —a{t)p'j (t), which is added in the right-hand 
side of Eq. HB2|) . Here a{t) is defined by 



a{t) = (B4) 

as a function of (t) and p'j (t) so that the total kinetic 

energy is constant in time: d\J2'j'=i |Pj(*)P /{2m)]/dt — 
0. Note that adding the thermostated term —a{t)p'j{t) 
we ignore Galilean invariance of the dynamics. There- 
fore the Sllod dynamics with the isokinetic thermostat 
is usually formulated in the center of mass frame and 

Ef=ip;w = o. 

Now, we introduce q'j^ il'jy) ^ the 1st component (2nd 
component) of the coordinate variable q^ of the j-ih par- 
ticle, and p'j,j. ijp'jy) as the 1st component (2nd compo- 
nent) of the variable of the j-ih particle. We define 
the operator L' by 



XT ( dU(q:(t)) 



iL'X'[T') 

N 

= E 



dX'{T') dH'{T') dX'{T') dH'{T') 



aq' 



9p' 



5q' 



N 



dX'{T') , dX'{T') , 



^3V 



dp'- 



(B5) 



for any function X'{T') of r'(= (p',q'))- Here H'{T') is 
introduced as 



^ I /|2 



(B6) 



As an important feature of the operator iV we have 



iL'H'iV) = -7VP,.y(r') 



(B7) 



using the function Pxy{T) defined by Eq. (|37|l . Using 
the operator iL' defined by Eq. IjBSp . the Sllod dynam- 
ics (jBl|l and (|B2|) is simply represented as dT'{t)/dt = 



iL'r'{t) for r'{t) = exp{iL'{t - to)}r' = (p'(i),q'(t)), 
which is equivalent to the equation 



dg'{T\t) 



dt 



= -zLV(r',t) 



(B8) 



for the distribution g'(T' , t) of F' at time t, because of the 
relation /dr'X(rV(r',t) = jdT'X{T'{t)) g'{T',to) 
for any function X(r') of F'. In other words, Eq. (|B8p is 
the Liouville equation corresponding to the Sllod equa- 
tion (|BlJ and (|B2)| . 

We define the distribution functions 5o(r') and g'(r,t) 

as 



5^(r')^S'-iexp{-/3i7'(r')} (B9) 

g'(r',t) = exp{-zL'(i-io)}5(,(r') (BIO) 

= S'expj-/? H'{T')+ ds 

9exp|-iL'(s-io)}i?'(r') 
ds 

= 5^,(r')exp|-/37V^ dsP;y(r',-s + 2to) 

(Bll) 

where P^y(r', t) is defined by by P^y{T\ t) = exp{iL'(t- 
io)}-Pxy(r') and S' is a normalization constant S' = 
/dr'exp{-/3i/'(r')}. Here we used the relation (|B7j) 
to derive Eq. IjBlip . Using the distribution f{T,t) we 
define the average (X(r'))j by 



(x(r'));^ drx{r')g'{T',t) 



(B12) 



for any function X(r') of T' . 

Now we discuss a relation between the above Sllod 
dynamics approach and the Hamiltonian dynamics ap- 
proach of the text of this paper. First, in the introduc- 
tion of the Sllod equation we notice correspondences of 
the spatial coordinate and the velocity as 



(mov) 



(B13) 
(B14) 



where the quantities on the left-hand side are for the 
Sllod dynamics approach, and the ones on the right-hand 
side are for the Hamiltonian approach of the text. Sec- 
ond, it is clear that Eqs. (|B9|) and (|B11|) correspond to 
Eqs. H42|l and H47|) . respectively, therefore we also have 
correspondences as 
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5o(r') 

iL' 

{■■■Yt 



/(r) 

■£(.me) 

/(r,i) 



(B15) 
(B16) 
(B17) 
(B18) 



These correspondences are not mathematical equiva- 
lences, but we can discuss some physical meanings in 
them. As an example, let's consider the correspon- 
dence (|B15|) more concretely. We introduce the new 



variables p' 



7eq', 



and 



p' = {p[,p2,--- ,p'n)- Here, T is the vector for the 
SUod equation, which corresponds to T for the approach 
of the text: T < — > T. Using this vector F we rewrite 
the distribution go{T') as 



N 




(B19) 



where we used the function form (r) of T given 

by Eq. (|^ . The Shod dynamics approach does not take 
into account the inertial force, as same as the approach 
using the distribution H13|l . and this is the reason why 

the second term J2j=i "^(7'?jy)^/2 in the square bracket 
of the right-hand side of Eq. (|B19|I . which make a dif- 
ference of g'o{r') from /(F) = S^^ exph/^TJ^™""' (F)], 
appears. However as far as we can neglect the effect of 
the inertial force, the distribution g'^iT') can correspond 
to the distribution /(F) in the shear flow system, namely 
we obtain the correspondence HB15|I . 

It should be noted that the vector F' used in the SUod 
equation does not correspond to the phase space vec- 
tor F, so the quantity Prcy{T') on the right-hand side of 
Eq. HB7|I does not generally correspond to the quantity 
Pxy{T) on the right-hand side of Eq. (jlTT) . Actually, 
noting the corresponding l|B13|l we notice 



N 



P-y(^') + Y,Y.1>'Ty ^ ^-^(r). (B20) 



However, if the two quantities q'jy and p'jy are decou- 
pled, namely {q'jyP'jy)'oc ~ WjyYocip'jyYoc^ and the average 
{p'jyYoa is zero, then we can have an approximate corre- 
spondence {P^ylT'))'^ < . {Pxy{T))oo- 

We can also show that using the distribution g'{T',t) 
based on the SUod equation, we can calculate the viscos- 
ity 77 as 



(P™(F'))' 
^ = _ lim ^ '^^^ 
7^0 7 



~ — lim 



JdT'P„/T')g',(T') 



lim/3V / ds f dT'P^y(T') 
xP^^(F',-s + 2io)5o(r') 



= /3V / dt {P,y{T,t)P^y{T) 

•I to 



(B21) 



(p',q') with where we used 



j dF'P,,(F')5o(r') = (B22) 

whose derivation is similar to that of Eq. (|52|l . noting 
the correspondences (|B15|) and (jB16|) . Therefore, using 
the SUod dynamics approach we obtain the same linear 
response formula for viscosity as Eq. 1)51(1 . 



APPENDIX C: SECOND LAW OF 

THERMODYNAMICS IN THE 
NON-EQUILIBRIUM CANONICAL 
DISTRIBUTION APPROACH 

In this appendix we give a derivation of the inequality 
(|^ satisfied at any time t (> to)- 
We start our derivation from the inequality 



a;lna:: — a:: + l>0 



(CI) 



satisfied by any positive real number x (> 0). The equal- 
ity in l|Cip is satisfied only when x = 1. Using the in- 
equality (|Cip in the case x = f{T, t)/ f{T) we have 



^^^^'\nim^^ + l>0, (C2) 



/(r) /(F 

which is equivalent to 



/(r) 



/(F, t) ln/(F, t) - /(F, t) In /(F) > /(F, t) - /(F) 



Now we note 



dTf{T,t) 



dF/(F) (= 1) 



(C3) 



(C4) 



dTfiTA) \nf(r,t) 
= J dFe-^^*'"^'(*-*°) [/(F) In /(F)] 

dF/(F)ln/(F). (C5) 
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By taking integral with respect to T on both sides of the 
inequahty (|n3|) . and by using Eqs. iglj, and (lUs)) 

we obtain 



I dTf{r,t)S{T)> I dTf{T)S{T). 



(C6) 



Using the equation S = /dr/(r)S'(r) = {S)to in the 
inequahty HC6|I . we obtain the inequahty lfS7|) . 



APPENDIX D: STABILITY CONDITION FOR 
THE SHEAR FLOW 

In this appendix we show the equivalence between the 
condition H66|l and the conditions (|67|) and Ht)8|) . We also 
give a derivation of Eq. . 

Noting that the energy is the function of S and 

Q by Eq. we have 



6T ^ 6 ^ ^ (55" + SQ (Dl) 



95 

Using Eq. (|DlJ we also have 



dSdQ 



dS 



dSdQ 



-6S5Q 



=2— W 

dQ 



^ dSdQ J dQ 



{6Tf 



Q2jj{ine) 



dSdQ 



dT 



{5Tf 



9^ 



{5Qf 
(D4) 



The inequality 1)66(1 must be satisfied by any infinitesimal 
deviations 5T and 5Q, so using Eq. (|D4p we obtain the 
conditions 1)67(1 and ((68(1 . 

Now, using the canonical distribution ((42(1 . we calcu- 
late the derivative of Q with respect to 7 at constant 
temperature T. 



07 = 



-SS ^ — SQ 



dSdQ 

Q2jj{me) 

dSdQ 



dQ 



g2jj{ine) 



g2jj{ine) _ 

ST - SQ 

dSdQ 



_2 SQ 
9Q 



dQ" 



g2jj{ine) 

ST 

dSdQ 



g2jj{ine) / g2jjit 



g2jj{me) 

dSOQ 



dQ 
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I^S- / dTQiT) 



exp{-/3 [H(™^)(r)-7Q(r)]} 



/ dT[QiT)Y 



exp{-/3 [i/(-^)(r)~7Q(r)]} 



dj 



dTQiT) 



exp{-/3 [i?(™^)(r)-7Q(r)]} 



(D5) 



where we used ^~^d:i/dj — j3Q. Therefore we obtain 
(D2) Eq. m- 



which leads to 



97 
9^ 



9Q 



g2jj{ine) / g2ff{ine) 



dS 



g2ff{ine) 

dSdQ 



Using Eq. ((Dip and ((D3(l we obtain 



APPENDIX E: RELATION BETWEEN THE TWO 
2 AVERAGES 

In this appendix we give a derivation of Eq. 1(71(1 . 

(D3) Using the expression 113) for the distribution f{T,t) 

used in the average {X{T))t for any function XiT), we 
have 
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^X + {{X-X))^ 
_ /•* d((X~X)) 



to 



ds 



X-p-fV f ds( [X(r) - X] P„j{T, -s + 2to: 

Jtn ^ 



to 

= X-p-fV [ ds 

J to 
nt 



X{T,s)~X 



X-PJV ds X{T, S)~X P,y{T) - P^y 

J to 



(El) 



where we used Eqs. (|5^ . {{X — X))t„ — and 
(iL(™'=))t = By taking the hmit < ^ oo in 

Eq. HEll) . we obtain Eq. (|7T|l . Concerning Eq. l|Eip one 
may notice 



X{T,s)-X P„j{T)-P,y =X{T,s)P,y{T 



(E2) 



because of Eq. (|52|l . so the integral function in the second 
tern of the right-hand side of Eq. IjElj) can be replaced 
by the right-hand side of Eq. (|E^ . 
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